using Pkg Pkg.activate(".") using DataFrames, CSV, StatsPlots, Downloads, DataFramesMeta, Dates, Statistics, Turing, Memoization, ReverseDiff, Serialization
if !Base.Filesystem.ispath("data/nics-firearm-background-checks.csv") Downloads.download("https://github.com/BuzzFeedNews/nics-firearm-background-checks/raw/master/data/nics-firearm-background-checks.csv","./data/nics-firearm-background-checks.csv") end if ~Base.Filesystem.ispath("data/fipscodes.csv") Downloads.download("https://www2.census.gov/geo/docs/reference/state.txt","./data/fipscodes.csv") end fipscodes = @chain CSV.read("data/fipscodes.csv",DataFrame) begin @subset(:STATE .< 60) end # manually downloaded CDC wonder gun homicide data using online form gunhomicidesnew = CSV.read("./data/wonder-gun-homicide-byyear.csv",DataFrame) gunhomicidesold = CSV.read("./data/cdc-wonder-firearm-homicide-1979-1998.csv",DataFrame) gunhomicides = [gunhomicidesnew ; gunhomicidesold] gunhomicides = @chain rename(gunhomicides,Dict("Crude Rate" => "crudedeathrate", "State Code" => "StateCode", "Year" => "year", "Year Code" => "YearCode")) begin @subset(.! ismissing.(:Deaths)) @transform(:crudedeathrate = map(x -> isnothing(x) ? missing : x, tryparse.(Float64,String.(:crudedeathrate))), :state = String31.(:State)) @subset(.! ismissing.(:crudedeathrate)) @orderby(:StateCode,:year) end gunchecks = CSV.read("./data/nics-firearm-background-checks.csv",DataFrame) startpops = @chain gunhomicides, @subset(:year .== 1999), @select(:State, :StartingPop = :Population,:StartDeathRate = :crudedeathrate) gunsalesest = @chain gunchecks, @orderby(:state,:month), @select(:state,:month,:salesest = 1.1 .* (:handgun .+ :long_gun) .+ 2.0 .* :multiple), groupby(:state), @transform(:cumsales = cumsum(:salesest),:year=year.(:month)), @subset(month.(:month) .== 6) joined = leftjoin(gunhomicides,gunsalesest, on=[:State => :state,:year => :year]; matchmissing=:notequal,makeunique=true) joined = @chain leftjoin(joined,startpops; on=[:state => :State],makeunique=true) begin @subset(.! ismissing.(:Deaths) .&& .! ismissing.(:StartingPop)) @subset(:StateCode .< 60) @transform(:cumgunrate = (:cumsales .+ 1.0 .* :StartingPop) ./ :Population ) end joined = @orderby(leftjoin(joined,fipscodes,on = :state => :STATE_NAME),:StateCode,:year) display(plot(joined.cumgunrate,joined.crudedeathrate; group=joined.state, legend=false, xlab="Gun Rate (guns/capita)", ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Gun Prevalence", alpha=.5)) display(plot(joined.year,joined.crudedeathrate; group=joined.state, legend=false, xlab="Year", ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Time", alpha=.5))
How has "Constitutional Carry" affected violence in states that pass those laws? Let's take a first look at the data, by plotting the firearms homicide rate through time, with a marker for the onset of Constitutional Carry as found in the Wikipedia article on Constitutional Carry.
concarry = leftjoin(fipscodes,CSV.read("./data/ConstCarryDates.csv",DataFrame),on = :STUSAB => :StateCode) concarry.date = map(x -> ismissing(x) ? missing : Date(div(x,10000),div(mod(x,10000),100),1),concarry.LawDate) concarry.year = map(x -> ismissing(x) ? missing : year(x),concarry.date) concarry = @subset(concarry,:STATE .<= 59) #ignore minor territories let p = [] for (st,stusab,statename,statens,lawdate,stdate,year) in eachrow(concarry) d = year sub = @subset(joined,:state .== statename) if(nrow(sub) > 0) pp = plot(sub.year,sub.crudedeathrate,ylim=(0,15),size=(500,500), title="$statename death rate vs time",xlab="Year",ylab="Gun Homicide/100k/yr",legend=false) if(!ismissing(d)) pp = plot!([(d,0),(d,15)]) end push!(p,pp) end end plot(p...; size = (2000,2000)) end
Looking at the raw data we can see that there are a number of states that have an "accelerating" trend of increased firearms homicide violence in the period from about 2015 onward. Some of these have converted to Constitutional Carry near the beginning of that period or during that period, with the trend beginning before the carry provision for several of the states. Also many states have passed their laws in the last few years and have no data post-law.
Let's look at these states grouped into geographic regions, because trends in regions may be more obvious. We will start with the following groupings:
stategroups = [["WA","OR","CA","NV","AZ"],["ID","MT","WY","UT","CO","NM","TX"],["ND","SD","NE","KS","OK","IA","MO"], ["MN","WI","IL","IN","MI","OH"],["AR","LA","MS","AL","GA","FL","SC"],["TN","NC","KY","WV","VA","MD","DE"], ["PA","NJ","NY","CT","RI","MA","VT","NH","ME"],["AK","HI","DC"]] stategroupdf = DataFrame(STUSAB = reduce(vcat,stategroups),stgroup = reduce(vcat,[[k for j in 1:length(stategroups[k])] for k in 1:length(stategroups)])) stategroupdf = leftjoin(fipscodes,stategroupdf; on = :STUSAB) @assert(length(unique(reduce(vcat,stategroups))) == 51) # 50 states plus DC for statelist in stategroups sub = joined[in.(joined.STUSAB, Ref(statelist)),:] ccdates = @chain @subset(concarry, in.(concarry.STUSAB,Ref(statelist))) begin @subset(.! ismissing.(:year)) end ccdatesenact = leftjoin(ccdates,sub, on = [:STUSAB => :STUSAB, :year => :year], makeunique = true) ccdatesenact = @subset(ccdatesenact, .! ismissing.(:crudedeathrate)) #display(ccdatesenact) #display(sub) p = plot(sub.year,sub.crudedeathrate,group=sub.STUSAB,ylim=(0,20),legend=:topleft,linewidth = 3, thickness_scaling=1, title="Firearm Homicide Rate vs Time", ylab="Homicides/100k population/yr") if nrow(ccdatesenact) > 0 p = scatter!(ccdatesenact.year,ccdatesenact.crudedeathrate,markersize=6,group = ccdatesenact.STUSAB) end display(p) end
Each state has its own overall level of firearm violence, but often trends in the same direction relative to its region even if its overall level is higher or lower. This suggests a model in which we measure each state using a dimensionless number relative to some baseline level. For example the average rate in the years 1999,2000,2001
baselines = @chain joined begin @subset(in.(:year,Ref((1999,2000,2001)))) @by(:STUSAB, :baseline = mean(:crudedeathrate)) end bljoined = @chain leftjoin(stategroupdf,leftjoin(baselines,joined,on = :STUSAB, makeunique=true),on=:STUSAB,makeunique=true) begin @subset(.! ismissing.(:crudedeathrate)) @orderby(:stgroup,:STUSAB,:year) end function getccd(gr) ccd = innerjoin(gr,@subset(concarry,.! ismissing.(:year)), on = [:STUSAB => :STUSAB, :year => :year], makeunique = true) ccd = @subset(ccd, .! ismissing.(:crudedeathrate)) ccd end function makerelplot(gr) ccd = getccd(gr) # display(ccd) p = plot(gr.year,gr.crudedeathrate ./ gr.baseline,group=gr.STUSAB,legend=:topleft,ylim=(0.0,2.5),xlab="Year",ylab="Relative Rate") if nrow(ccd) == 1 scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline) elseif nrow(ccd) > 1 scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline,group = ccd.STUSAB) end p end
makerelplot (generic function with 1 method)
Each state plotted in a group with other geographically similar states. Firearm homicide rate relative to the recorded rate in the year 1999.
plots = [makerelplot(gr) for gr in groupby(bljoined,:stgroup)] plot(plots... ; size=(1000,1000),linewidth=2,xlab="Year",ylab="Relative Rate")
We can try to estimate the effects of these laws by building a model of the overall process. Because we can't resurrect people from the dead, the homicide rate can never be negative. It makes sense then to model the homicide rate on a logarithmic scale.
We estimate the behavior of the states as an overall level, plus a shape that is common to the region. Individual states are then allowed a small perturbation to the regional shape.
The effect of constitutional carry law would then be estimated as the log firearm homicide rate which exceeds the state level counterfactual estimate in the years after the passage of the law.
We impose an informal smoothness requirement on the counterfactual estimates by using a compact radial basis function expansion with one center every 4 years, and a maximum radius of 7 years so that there is overlap between adjacent centers.
function bumpfun(x,c,scale) stdx = (x-c)/scale if stdx < -1.0 || stdx > 1.0 0.0 else exp(1.0-1.0/(1.0-stdx^2)) # goes to zero at -1 and 1, and 1 at x=0 end end function timeser(x,coefs,centers,scale) f = 0.0 for (a,c) in Iterators.zip(coefs,centers) f += a*bumpfun(x,c,scale) end return f end function laweffect(yr,rate,start) if yr >= start 1.0-exp(-rate*(yr-start)) else 0.0 end end statecoefpert = log(1.05)/2.0 include("model.jl") ## load the model, this is a separate file so we can save the output unless it changes modeljoined = @chain leftjoin(joined,stategroupdf; on = :STUSAB, makeunique=true) begin @subset(.! ismissing.(:crudedeathrate)) end lawdates = @chain leftjoin(DataFrame(STATE=1:55),concarry; on=:STATE) begin @subset(:STATE .< 60) @rtransform(:year = ismissing(:year) ? 3000 : :year) @orderby(:STATE) end centers = [i for i in (minimum(modeljoined.year)-4):4:(maximum(modeljoined.year)+4)] width = 7.0 modl = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup), centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),lawdates.year,statecoefpert) setadbackend(:reversediff) Turing.setrdcache(true) savedfile = "./saved/samples.dat" global s = [] if stat(savedfile).mtime > stat("./model.jl").mtime global s = deserialize(savedfile) else global s = sample(modl,NUTS(500,0.8),MCMCThreads(),500,3) serialize(savedfile,s) end
Having sampled the model, we can then compute summary graphs of the results. We will show the posterior density for the coefficient of the magnitude of the law effect for each state separately.
Also we will plot the actual log(homicide rate) together with the state specific counterfactual estimates.
lawco = group(s,"statelawcoef") statenames = Dict() for (st,stusab) in Iterators.zip(lawdates.STATE,lawdates.STUSAB) statenames[st] = stusab end function plotlawts(s,lawdates,modeljoined) let pl = [] lawco = group(s,"statelawcoef") for i in @select(@subset(lawdates,:year .< 2022),:STATE).STATE den = density(s[:,Symbol("meanlawcoef"),:] .+ lawco[:,Symbol("statelawcoef[$i]"),:], title = "State $i = $(statenames[i])",xlim=(-log(2),log(2))) plot!([(0.0,0.0),(0.0,5.0)], color="red",ylim=(0.0,5.0)) push!(pl,den) end println("State law coefficient values") display(plot(pl...; size = (1000,1000))) display(density(s[:,:lawrate,:], title= "Law effect onset rate (1/yr)")) modeljoined.logdrate = log.(modeljoined.crudedeathrate) regions = Dict() for r in eachrow(stategroupdf) regions[r.STATE] = (group=r.stgroup,code=r.STUSAB) end pl = [] samps = sample(1:500,10) for st in unique(modeljoined.STATE) sub = @subset(modeljoined, :STATE .== st) region = regions[st].group p = plot(sub.year,sub.logdrate,linewidth=3,title="$(st) = $(regions[st].code)",ylim=(-0.5,2.75)) push!(pl,p) for samp in samps statelawcoef = s[samp,Symbol("statelawcoef[$st]"),1] statecoefs = [s[samp,Symbol("statecoefs[$st][$i]"),1] for i in 1:length(centers)] statebase = s[samp,Symbol("statebase[$st]"),1] regioncoefs = [s[samp,Symbol("regioncoefs[$region][$i]"),1] for i in 1:length(centers)] lawrate = s[samp,:lawrate,1] startlaw = lawdates[lawdates.STATE .== st,:].year[1] startlaw = ismissing(startlaw) ? 3000 : startlaw years = 1979:2020 pred = [statebase + timeser(yr,regioncoefs .+ (statecoefpert .* statecoefs),centers,width) for yr in years] plot!(years,pred; color="orange",alpha=0.5) if startlaw < 3000 plot!([(startlaw,-10.0),(startlaw,10.0)],color="red") end regpred = [statebase + timeser(yr,regioncoefs ,centers,width) for yr in years] end end println("Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)") display(plot(pl...; size =(2000,3000), legend=false)) end density(s[:,:meanlawcoef,:],title="Mean law coefficient") end plotlawts(s,lawdates,modeljoined)
State law coefficient values Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)
One plausible explanation for the relative lack of clear measurable effect from permitless carry laws is that permitless actually may have had very little effect on how many people carry firearms in public. Many of these states were shall-issue and many people who wanted to carry may have had CCW permits issued before the permitless law passed. Therefore the change in number of people carrying may have been minimal.
Using the date on which states converted to shall-issue permits, with a similar analysis, we can determine whether more of an effect was visible after that legal change. Some states converted from no-issue to shall issue, while others converted from may-issue to shall-issue.
# dates collected from https://en.wikipedia.org/wiki/History_of_concealed_carry_in_the_United_States shallissdates = CSV.read("data/ShallIssueDates.csv",DataFrame) shallissue = leftjoin(fipscodes,shallissdates; on = :STUSAB) silawdates = @chain leftjoin(DataFrame(STATE=1:55),shallissue; on=:STATE) begin @subset(:STATE .< 60) @rtransform(:year = ismissing(:ShallIssueDate) ? 3000 : :ShallIssueDate) @orderby(:STATE) end modl2 = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup), centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),silawdates.year,statecoefpert) savedfile = "./saved/shallisssamples.dat" global sisam = [] if stat(savedfile).mtime > stat("./model.jl").mtime global sisam = deserialize(savedfile) else global sisam = sample(modl2,NUTS(500,0.8),MCMCThreads(),500,3) serialize(savedfile,sisam) end
plotlawts(sisam, silawdates,modeljoined)
State law coefficient values Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)
Many states that have passed constitutional carry laws have done so recently enough that there is no data on homicides available from CDC Wonder yet.
For those which passed the law long enough ago to have some post-law data, still many of them were in the last few years leaving only a few years of post-law data. The last few years is a time period where homicides have increased nearly nationwide and it is challenging to estimate the effect of a law when it is imposed on top of a nonlinearly changing trend. Estimates of a law effect created by comparing the actual data to counterfactuals estimated as small perturbations to the regional trend lead to estimates of the average law coefficient that has uncertain sign, with a most likely value of about 0.05 or so, meaning that gun homicides may increase by a factor of around 1.05, however the ability to resolve this number is so poor that the number could be anything from about -0.1 to +0.25 (multiplying the crime rate by anywhere from 0.9 to 1.28.
Evidence for a consistent average net positive or negative effect of passing constitutional carry laws simply doesn't exist in the data so far.
In states such as AK, AZ, and AR which passed their laws long enough ago that we might expect to be able to estimate the effect, the average effect appears to be close to 0.0 with uncertainty relatively symmetrically distributed to both sides. The strongest estimates we have are near zero.
KY, and MS appear to be outliers, but even in these cases it's ambiguous whether the law resulted in an increase or a decrease, even though increase appears to be the more probable direction.
The evidence simply does not show a consistent effect increasing or decreasing firearm homicides after liberalization of concealed carry laws.